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►  The  phase  shifts  and  time  scales  of  the  major  dynamic  processes  are  determined. 
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In  this  paper,  the  dynamic  response  of  a  molten  carbonate  fuel  cell  is  studied  by  means  of  computational 
fluid  dynamics.  First,  a  rigorous  three-dimensional,  transient  and  non-isothermal  model  is  developed 
through  a  comprehensive  inclusion  of  various  transport  phenomena.  Next,  a  sinusoidal  impedance 
approach  is  used  to  examine  the  dynamic  response  of  the  unit  cell  to  inlet  perturbations  at  impedance 
frequencies  of  1,  0.1,  0.01  and  0.001  Hz.  This  analysis  is  further  used  to  determine  the  phase  shifts  and 
time  scales  of  the  major  dynamic  processes  within  the  fuel  cell.  The  load-following  capability  of  the  unit 
cell  is  studied  by  examining  the  dynamic  responses  of  the  average  current  density,  electrochemical  re¬ 
actions  rates,  heat  and  mass  transfer,  mass  fractions  and  temperature.  The  results  show  that  the  elec¬ 
trochemical  reactions  and  the  charge  transport  process  adjust  within  a  millisecond.  The  mass  transport 
process  shows  a  comparatively  larger  time  scale  that  is  about  1  s.  The  heat  transport  process  is  the 
slowest  process  in  the  cell  and  takes  about  an  hour  to  reach  its  steady-state.  The  developed  model  along 
with  the  dynamic  analysis  results  can  be  used  to  design  process  control  strategies. 

©  2013  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Molten  carbonate  fuel  cells  have  been  intensively  studied  over 
the  past  few  decades,  and  many  technical  barriers  have  been 
overcome.  However,  cell  corrosion  and  lifetime  are  still  considered 
the  greatest  obstacles  towards  commercialization. 

From  an  application  perspective,  dynamic  operation  has  a  signif¬ 
icant  effect  on  the  cell  life  cycle,  and  hence  economic  matters.  The 
MCFC  behaviour  in  dynamic  situations  is  still  not  fully  understood. 
Specifically,  when  the  cell  undergoes  voltage/load  variation  or  fluc¬ 
tuation,  predicting  the  fuel  cell’s  dynamic  performance  becomes 
crucial.  The  transient  situation  at  MCFC  start-up  is  another  example 
for  varying  conditions.  Having  an  improved  understanding  of  the 
system  behaviour  at  dynamic  phase  helps  to  design  a  more  robust 
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control  system  in  order  to  minimize  fuel  consumption  and  maximize 
operating  lifetime.  In  fact,  in  order  to  choose  an  optimal  control 
system  and  operation  parameters,  an  appropriate  dynamic  model  of 
MCFC  consistent  with  other  components  of  the  hybrid  system  is 
ideal.  This  can  be  achieved  by  determination  of  time  scales  of  the 
major  dynamic  processes  in  the  MCFC,  and  their  analysis  and  com¬ 
parison  with  time  scales  of  other  devices  in  the  hybrid  system.  The 
developed  dynamic  model  can  be  further  employed  to  examine 
quantitative  responses  of  the  system  to  different  inlet  perturbations. 
In  addition,  phase  shifts  between  the  dynamic  responses  of  various 
system  characteristics  and  the  voltage  perturbation  will  be  identified. 

Over  the  past  three  decades,  several  studies  have  been  performed 
in  order  to  develop  comprehensive  mathematical  models  repre¬ 
senting  the  various  processes  which  occur  in  molten  carbonate  fuel 
cell  operation  [1-13].  Some  authors  have  implemented  these 
models  to  investigate  steady-state  performance  characteristics  of 
MCFCs  at  cell  and  stack  levels  [14-20].  Having  said  that,  imple¬ 
mentation  of  these  models  in  transient  simulations  is  often  subject 
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Nomenclature 

e 

porosity 

6 

electrolyte  filling  degree 

Ai 

electrode  active  surface  area  (m2  nrr3) 

V 

overpotential  (V) 

cp 

specific  heat  (J  kg-1  K-1) 

F 

dynamic  viscosity  (kg  nrr1  s_1) 

D 

mass  diffusivity  (m2  s-1) 

V 

species  stoichiometric  coefficient  of  the  reaction 

£eq 

equilibrium  electric  potential  (V) 

P 

density  (kg  nr3) 

F 

Faraday’s  constant,  96,485  (C  mol-1) 

(p 

electric  potential  (V) 

f 

frequency  (Hz) 

a 

electric  conductivity  of  the  solid  phase  (S  nr1) 

io 

exchange  current  density  (A  m-2) 

K 

electric  conductivity  of  the  electrolyte  phase  (S  m_1) 

i° 

*0 

reference  exchange  current  density  (A  nr2) 

j 

local  current  density  (A  m  2) 

Subscripts  and  superscripts 

R 

volumetric  current  density  (A  nrr3) 

a 

anode 

k 

thermal  conductivity  (W  nr1  I<-1) 

age 

anode  gas  channel 

I< 

permeability  (m2) 

c 

cathode 

M 

molecular  weight  (g  mol-1) 

ege 

cathode  gas  channel 

n 

number  of  electrons 

e 

electrolyte  phase 

P 

static  pressure  (Pa) 

g 

gas  phase 

R 

universal  gas  constant  (8.314  J  mol-1  I<-1) 

in 

inlet 

S 

molar  entropy  (J  mol-1  I<-1) 

i 

ith  component 

S 

source  terms 

j 

jth  species 

t 

time  (s) 

m 

mass  equation 

T 

temperature  (I<) 

out 

outlet 

V 

gas  velocity  (m  s-1) 

s 

solid  phase 

X 

molar  fraction 

T 

energy  equation 

y 

mass  fraction 

u 

momentum  equation 

(j)e 

electronic  charge  equation 

Greek  letters 

(j)c 

carbonate  ion  charge  equation 

a 

transfer  coefficient 

eff 

effective 

T 

tortuosity 

ref 

reference  state 

to  several  simplifications.  In  other  words,  model  reduction  methods 
are  widely  used  to  establish  a  simple  simulation  model  according  to 
various  operating  conditions.  For  instance,  Hao  et  al.  [21  ]  developed 
a  simplified  dynamic  model  for  a  cross-flow  molten  carbonate  fuel 
cell  and  solved  the  model  using  VC++.  The  model  was  two- 
dimensional  with  uniform  voltage  distribution.  This  model  was 
developed  based  on  identical  heat  characteristics  for  anode,  cathode 
and  electrolyte  although  they  typically  differ  by  an  order  of  magni¬ 
tude.  Another  two-dimensional  model  was  presented  by  Fermegila 
et  al.  [22]  to  study  the  effects  of  step  change  and  linear  ramp  for 
various  inlet  conditions.  The  authors  neglected  the  enthalpy  trans¬ 
port  between  the  electrode  pores  and  gas  phase,  and  also  assumed 
an  equally  distributed  current  pattern  in  both  electrodes  (which  is 
not  always  appropriate).  Likewise,  there  exist  few  studies  at  the 
stack  level.  For  instance,  Lukas  et  al.  [23]  employed  a  lumped- 
parameter  formulation  of  the  first  principle  equations  for  a  fuel 
cell  stack  with  a  simplifying  assumption  indicating  that  the  solid 
mass  temperature  is  equal  to  the  exit  stream  temperature.  In  fact, 
this  model  does  not  provide  sufficient  information  about  the  tem¬ 
perature  profile  in  the  electrodes  and  the  electrolyte  which  is  critical 
for  MCFC  lifetime  and  components  degradation. 

It  was  found  that  a  three-dimensional  dynamic  analysis  is  un¬ 
available  in  the  literature.  He  and  Chen  [24]  developed  a  three- 
dimensional  transient  stack  model  using  a  commercial  package  to 
demonstrate  the  heat  transport  at  the  stack  level.  However,  the 
processes  of  gas  transport  and  chemical  reactions  were  incorpo¬ 
rated  only  at  the  cell  level. 

The  dynamic  performance  of  fuel  cell  hybrid  systems  has  also 
been  investigated  by  some  researchers.  As  an  example,  an  inves¬ 
tigation  of  the  control  performance  of  an  internal  reforming  MCFC 
system  was  performed  by  Lukas  et  al.  [23,25,26].  The  authors  used  the 
actual  data  from  a  2-MW  demonstration  project.  However,  the  focus 
of  this  study  was  primarily  on  the  molten  carbonate  fuel  cell  and  the 


hybrid  system  was  not  considered.  Zhang  et  al.  [27]  presented  a  dy¬ 
namic  model  for  a  hybrid  fuel  cell-gas  turbine  system  with  some 
control  loops  applied  to  the  system.  A  distributed  power  generation 
system  was  studied  by  Grillo  et  al.  [28].  This  hybrid  system  is  based  on 
pressurization  and  heat  recovering  of  a  100  kW  molten  carbonate 
fuel  cell.  Au  et  al.  [29]  investigated  the  optimization  of  MCFC  oper¬ 
ating  temperatures  by  presenting  a  case  study  in  which  the  efficiency 
of  a  CHP  plant  was  analyzed.  Some  authors  [e.g.,  Ref.  [30]]  have  also 
used  energy  and  exergy  analyses  to  evaluate  various  system  effi¬ 
ciencies  for  the  integrated  power  generation  systems. 

A  neural  network  structure  approach  to  the  dynamic  perfor¬ 
mance  modelling  was  implemented  by  Shen  et  al.  [31]  to  model 
transient  cell  behaviour  of  a  fuel  cell.  However,  this  model  is  not 
appropriate  for  other  high-temperature  fuel  cell  systems. 

Recently,  a  transient  mathematical  model  for  a  single  counter¬ 
flow  MCFC  with  an  internal  reformer  was  developed  by  Heide- 
brecht  and  Sundmacher  [32].  Simplifications  like  plug  flow  and 
constant  pressure  in  the  gas  phase  along  with  a  lumped  solid  phase 
for  energy  balance  were  used.  Later  on,  this  group  [33]  presented 
a  more  detailed  study  based  on  a  dimensionless  mathematical 
model  of  a  single  cross-flow  MCFC  with  spatially  distributed  sim¬ 
ulation  results  for  steady-state  and  dynamic  scenarios.  This  model 
can  be  applied  to  any  other  high-temperature  fuel  cell  such  as  SOFC 
but  only  in  two-dimensional  simulation  mode. 

In  brief,  a  considerable  number  of  present  studies  have  been  car¬ 
ried  out  by  either  performing  the  simulation  for  individual  MCFC 
components  or  by  reducing  the  dimension  to  0D  (i.e.  black  box),  1 D  or 
2D.  In  addition,  many  approaches  have  employed  a  uniform  distri¬ 
bution  of  field  variables,  which  significantly  influences  the  accuracy  of 
the  dynamic  results.  Recently,  a  comprehensive  three-dimensional 
model  was  developed  in  Ref.  [34]  without  these  simplifying  as¬ 
sumptions.  However,  the  focus  was  mainly  on  MCFC  performance 
analysis. 
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Table  1 

The  geometric  parameters  of  the  simulated  MCFC. 


Parameter 

Value 

Cell  length,  (mm) 

100 

Cell  width,  (mm) 

100 

Anode  gas  channel  height,  (mm) 

2.0 

Anode  gas  channel  width,  (mm) 

2.0 

Cathode  gas  channel  height,  (mm) 

2.0 

Cathode  gas  channel  width,  (mm) 

2.0 

Anode  thickness,  (mm) 

0.7 

Anode  width,  (mm) 

4.0 

Cathode  thickness,  (mm) 

0.6 

Cathode  width,  (mm) 

4.0 

Electrolyte  thickness,  (mm) 

1.0 

Porosity  of  anode,  ea  [19] 

0.52 

Porosity  of  cathode,  ec  [19] 

0.62 

In  this  paper,  we  extend  our  previous  work  [35]  to  introduce 
a  rigorous  three-dimensional  transient  model  considering  a  broad 
incorporation  of  various  transport  phenomena.  The  main  objective 
of  this  research  is  to  assess  the  non-linear  dynamic  response  of  the 
unit  cell.  In  contrast  to  the  former  studies  that  have  used  ‘step 
change’  perturbation,  in  this  study,  a  sinusoidal  impedance 
approach  for  a  variety  of  frequencies  is  integrated.  In  fact,  there  is 
no  study  on  MCFC  dynamic  analysis  based  on  a  sinusoidal  impe¬ 
dance  approach  in  the  open  literature.  Based  on  this  analysis,  the 
dynamic  responses  of  various  field  variables  and  quantities  are 
measured.  Furthermore,  the  time  scales  of  various  transport  phe¬ 
nomena,  including  mass,  heat  and  electric  charge  transport  pro¬ 
cesses,  are  determined.  In  addition,  the  phase  shifts  for  the  non¬ 
linear  dynamic  response  of  field  variables  and  quantities  such  as 
temperature,  mass  fraction,  current  density,  etc.  are  identified. 

2.  Modelling 

A  planar  MCFC  is  a  rectangular  assembly  of  an  electrolyte 
sandwiched  between  an  anode  and  a  cathode.  The  fuel  gas  and 
oxidant  gas  are  introduced  through  anode  gas  channel  (AGC)  and 
cathode  gas  channel  (CGC),  respectively.  In  this  arrangement,  the 
fuel  gas  and  oxidant  flow  into  the  porous  anode  and  cathode  where 
the  hydrogen  oxidation  reaction  (Equation  (1))  and  oxygen  reduc¬ 
tion  reaction  (Equation  (2))  occur,  respectively: 

H2  +  COJ  <-►  H20  +  C02  +  2e~  ( 1 ) 

io2  +  C02  +  2e-~COf  (2) 

The  geometric  parameters  of  the  simulated  MCFC  are  summa¬ 
rized  in  Table  1. 


2.1.  Governing  equations 

A  comprehensive  MCFC  model  needs  to  consider  the  transport 
of  multi-component  gas  species  in  gaseous  and  liquid  phases, 
electrochemical  and  chemical  reaction  kinetics,  heat  generation, 
heat  transfer,  transport  of  electrons  and  carbonate  ions,  and 
porous  electrode  effects.  These  processes  occur  in  void  volumes, 
liquid  phase,  solid  phase  and  at  triple-phase  boundaries  (TPB). 
The  next  section  provides  a  brief  explanation  of  the  governing 
equations.  The  reader  is  referred  to  Ramandi  et  al.  [34]  for  more 
details. 

■  Transport  of  gas  species 

In  the  gas  phase,  the  governing  equations  include  mass,  mo¬ 
mentum  and  species  conservation: 

—  (£effpg)  +  (pg  ug)  =  Sm  (3) 

at(gS +V'  =  _VPg  +  V’P)+Su 

lt  (eefV!j  +  V •  (  -  PgDfmWi)  +  V •  (pg tg Yi)  =  Si  (5) 

where  pg  and  T tg  are  the  superficial  values  of  the  gas  mixture 
density  and  velocity,  respectively.  The  gas  mixture  density  (kg  m-3) 
is  calculated  based  on  the  ideal  gas  law  [36]: 

%  -  ^Eil)  '  <6> 

where  Pg  is  the  superficial  value  of  the  gas  pressure  (Pa),  T  the 
temperature  (I<),  R  the  universal  gas  constant  (J  kmol-1  I<-1).  In 
addition,  Y*  and  Mi  are  mass  fraction  and  molecular  weight 
(kg  kmol-1)  of  species  z,  respectively.  The  effective  porosity  (eeff)  is 
defined  using  the  filling  degree  of  electrolyte  and  implemented  in 
the  governing  equation  as  well  as  in  the  constitutive  laws: 

eeff  =  e(l  -  6)  (7) 

In  addition,  Sm  (kg  m-3  s-1),  Su  (Pa  m-3)  and  S*  (kg  m~3  s-1)  are 
the  mass,  momentum  and  species  source  terms,  respectively,  which 
have  different  values  with  respect  to  various  cell  sub-domains,  as 
given  in  Table  2.  Also,  Df^  represents  the  effective  diffusivity  of 
species  i  in  the  gas  mixture  and  is  evaluated  as 


Table  2 

Source  terms  in  various  conservation  equations. 


Sm 

S 

Su 

% 

St 

Gas  Channels 

0 

0 

0 

0 

0 

0 

Anode 

Sh2  +Sh2o  +  Sco2 

= 

-|M  * 

/(eff  U  g 

Ra 

Ra 

i  +  i+lvM*l^(TASa) 

Sh2o  = 

|mH!0 

Electrolyte 

0 

C 

II 

o 

(yf 

|mC02 

0 

0 

0 

Jl 

Cathode 

So2  +  Sco2 

i 

II 

H  T7 
j(e  ff  g 

Rc 

Rc 

K 

f+J+Mcl  +  ll'TASc) 

1 

O 

II 
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D! 


,eff 

i,  m 


l-Xj 

E  (Xj/Df ) 

ij^i  v 


(8) 


where  D|ff  is  the  effective  binary  mass  diffusion  coefficient  of 
species  i  in  species  j  as  calculated  by 


T  Pref  eeff 
neff  -  D  _  _ _ — 


(9) 


and  Dy  is  the  bulk  binary  diffusivity  at  the  reference  temperature 
(7*ef)  and  reference  pressure  (Pref).  Also,  t  is  the  tortuosity  of  the 
porous  material. 


■  Transport  of  electric  charge 


Fa  —  Av,a  '*0,a 


Xh2 


X, 


H2,in 


0.5 


exp 


«aaf 

RT 


’Ja  - 


-0.5 

Xh2  \  /  Xh2o 


X, 


H2,in 


Xco2 


X, 


C02,in 


exp(  “^’Ja 


X]-[20,in 

(16) 


Also,  the  superoxide  mechanism  [3]  is  employed  to  describe  the 
cathodic  reaction  rate: 


Rc  —  AV)c'io,c 


X, 


co2 


X, 


C02,in 

Xo2 


Xo2,in 


exp 

0.75 


a  acP, 

RT 


Vc  - 


Xt 


co2 


X, 


C02,in 


-0.5 


accF 


exp(  “"j^c 


(17) 


The  governing  equation  for  the  electronic  and  ionic  charge 
transport  in  MCFCs  can  be  derived  by  means  of  Ohm’s  law  as 
follows: 


where  Av,  F  and  a  are  the  electrode  active  surface  area,  Faraday 
constant  and  transfer  coefficient,  respectively.  i'o  is  the  exchange 
current  density  [19]  as  given  by 


on 

II 

-e? 

> 

fa 

CD 

b 

(10) 

AeffV4>e)  = 

(11) 

<0,a  =  i8.a(XH2,in)a25(XH20.in)°'25(Xco2.in)a25 

(18) 

l0,c  =  ^o,c(^02,in)°'625(Xco2,in)  °'75 

(19) 

The  overpotential,  7],  is  defined  as 

where  (7eff  and  /ceff  are  the  effective  electric  conductivity  of  the  solid 
material  and  liquid  electrolyte,  respectively.  These  two  parameters 
are  estimated  based  on  the  Bruggemann  correlation  [19]  as  follows: 


Va  =  0s  -  0e  (20) 

Vc  =  0s  “  0e  “  ^eq  (21 ) 


<r(l  -  e) 

(12) 

k(s8)15 

(13) 

The  source  terms  on  the  right  hand  side  of  Equations  (10)  and 
(11)  are  illustrated  in  Table  2. 


Moreover,  Feq  is  the  potential  difference  between  solid  and 
electrolyte  phase  potentials  in  equilibrium,  i.e.  when  no  current  is 
generated,  and  is  defined  using  the  Nernst  equation  [34]: 


Feq  —  Fq  + 


RT,  (PH2.aPc02.cPo)c\ 

nF  \  PC02,aPH20,a  ) 


(22) 


■  Transport  of  energy 

The  energy  equation  applying  to  each  sub-domain  of  the  fuel 
cell,  can  be  written  as  [37]: 

£  (^Cp)kTN)+V-(-/ceffVT)+V-(Tfg£pgCpT)=ST  (14) 
\k= g,s,e  / 

where  keff  is  the  effective  thermal  conductivity  determined  by 

kefr  =  (1  -  e)ks  +  e(l  -  6)kg  +  eeke  (15) 

Flere,  /<s,  kg  and  ke  are  the  thermal  conductivity  of  the  solid 
material,  gas  mixture  and  liquid  electrolyte,  respectively. 

The  heat  generation  or  consumption  is  represented  by  the 
source  term,  Sj.  Three  kinds  of  heat  sources  are  considered  in  the 
presented  model,  namely  the  reversible  heat  release  during  the 
electrochemical  reaction,  irreversible  or  activation  heat  generation 
and  ohmic  heating.  Table  2  summarizes  the  values  of  energy  source 
terms  in  each  sub-domain. 

■  Electrochemical  reaction  rates 

All  of  the  governing  equations,  described  in  the  above  sections, 
are  strongly  coupled  and  dependent  on  the  electrochemical  reac¬ 
tion  rates,  according  to  Table  2.  For  the  anodic  reaction  rate,  the 
correlation  proposed  by  Ang  and  Sammels  [38]  is  incorporated  as 


£0  =  1.2723  -  2.7645  x  1(T4T  (23) 

Table  3 

Physical  and  thermal  properties  of  various  materials. 

Parameter  Value 

Thermal  conductivity  of  anode,  (W  m  1  I<-1)  [15,43]  78 

Thermal  conductivity  of  cathode,  (W  m  1  K-1)  [15,43]  0.9 

Thermal  conductivity  of  electrolyte,  (W  m  1  K-1)  [15]  2.0 

Specific  heat  of  anode,  (J  kg-1  I<-1)  [44]  444 

Specific  heat  of  cathode,  (J  kg-1  K-1)  [44]  4435 

Specific  heat  of  electrolyte,  (J  kg”1  I<-1)  [44]  4000 

Density  of  anode,  (kg  m”3)  [44]  8220 

Density  of  cathode,  (kg  m-3)  [44]  6794 

Density  of  electrolyte,  (kg  m-3)  [44]  1914 

Electric  conductivity  of  anode,  (S  m_1)  [9,45]  1300 

Electric  conductivity  of  cathode,  (S  m”1)  [9,45]  1300 

Free  electrolyte  conductivity:  pre-exponential  factor,  (S  m_1)  [2]  3637 

Free  electrolyte  conductivity:  Apparent  activation  energy,  (K_1)  [2]  3016 

Hydrogen  diffusivity  in  carbon-dioxide,  (m2  s_1)  [44,46]  5.5E-5 

Hydrogen  diffusivity  in  water  vapour,  (m2  s_1)  [10,44]  9.15E-5 

Oxygen  diffusivity  in  carbon-dioxide,  (m2  s_1)  [9,46]  1.4E-5 

Oxygen  diffusivity  in  nitrogen,  (m2  s_1)  [44,46]  1.8E-5 

Carbon-dioxide  diffusivity  in  water  vapour,  (m2  s_1)  [44,46]  1.62E-5 

Carbon-dioxide  diffusivity  in  nitrogen,  (m2  s_1)  [44,46]  1.6E-5 

Hydrogen  diffusivity  in  liquid  electrolyte,  (m2  s_1)  [47]  IE-7 

Oxygen  diffusivity  in  liquid  electrolyte,  (m2  s_1)  [36,47]  3E-7 

Carbon-dioxide  diffusivity  in  liquid  electrolyte,  (m2  s_1)  [36,47]  IE-7 
Water  vapour  diffusivity  in  liquid  electrolyte,  (m2  s_1)  IE-7 

Nitrogen  diffusivity  in  liquid  electrolyte,  (m2  s_1)  IE-7 

Standard  entropy  change  of  anode,  (J  mol-1  K-1)  [48]  54.56 

Standard  entropy  change  of  cathode,  (J  mol-1  K-1)  [48]  -216.2 

Standard  entropy  change  of  generation  reaction,  (J  mol”1  K”1)  [48]  -161.64 
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Fig.  1.  Boundary  conditions  for  the  MCFC  model:  (a)  down  the  channel  view,  (b)  cross-sectional  view. 


The  physical  and  thermal  properties  that  are  employed  in  the  ■  At  AGC  inlet  (/agc)  and  CGC  inlet  (/cgc) 
developed  model,  are  listed  in  Table  3. 

The  gas  species  mass  fraction,  total  mass  flux  and  temperature 
2.2.  Boundary  and  initial  conditions  of  the  entering  gas  flow  are  specified.  Furthermore,  the  normal 

fluxes  of  all  other  variables  are  set  to  zero. 

(i)  Boundary  conditions 


A  complete  set  of  boundary  conditions  is  required  for  the  gov¬ 
erning  equations  formulated  in  the  previous  section.  It  should  be 
specified  involving  the  transport  of  gas  species,  energy,  electrons 
and  carbonate  ions.  A  demonstration  of  the  various  internal  and 
external  boundaries  is  presented  in  Fig.  1.  The  boundary  conditions 
are  required  only  at  the  external  surfaces  of  the  computational 
domain  due  to  the  implemented  single-domain  formulation 


Yf  =  Specified 
m-Yt  =  Specified 
T  =  Specified  > 

^  =  0 
67? 


(24) 


Fig.  2.  Illustration  of  the  modified  solution  procedure  based  on  the  developed  in-house  code,  implemented  in  ANSYS  FLUENT. 
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Fig.  3.  Model  parameter  verification:  using  the  experimental  study  by  Lee  et  al.  [39]. 


0  0.25  0.5  0.75  1 

Non-Dimensional  Time 

Fig.  5.  The  applied  sinusoidal  voltage  perturbation  with  amplitude  of  0.1  V  under 
different  impedance  frequencies  (10.0-0.001  Hz)  during  one  sinusoidal  cycle.  The 
horizontal  axis  represents  dimensionless  time. 


(26) 


■  At  anode  lands  (Wa>L) 


©  represents  any  remaining  variable  which  is  not  explicitly 
specified. 

■  At  AGC  outlet  (Oagc)  and  CGC  outlet  (Ocgc) 

Since  the  gas  channels  aspect  ratio  is  very  large,  the  flow  is 
assumed  to  be  fully  developed.  This  means  that  none  of  the  vari¬ 
ables  and  respective  fluxes  varies  in  the  normal  direction.  In 
addition,  the  gas  pressure  is  specified. 


P  =  Specified 

a© 


On 


&  0 


(25) 


The  electronic  potential  is  set  to  zero.  The  flux  of  ionic  charge 
and  all  remaining  variables  are  set  to  zero. 


(27) 


■  At  cathode  lands  (Wc,l) 

The  electronic  potential  is  set  to  be  the  cell  voltage.  The  flux  of 
ionic  charge  and  all  remaining  variables  are  set  to  zero. 


(28) 


■  At  AGC  walls  (Wage)  and  CGC  walls  ( Wcgc) 

A  no-slip  boundary  condition  is  applied  to  these  walls  along 
with  a  zero-flux  boundary  condition  for  all  other  variables. 


■  At  no-flux  boundaries  (WL  and  W^) 

The  zero-flux  condition  is  assumed  for  all  variables. 

{H  -  °}  ,29) 

(ii)  Initial  conditions 

For  the  case  of  sinusoidal  voltage  perturbation,  the  initial  con¬ 
ditions  are  considered  to  be  the  steady-state  simulation  results 
with  the  same  model  input  parameters. 


Fig.  4.  Model  validation:  using  the  experimental  and  numerical  study  by  Brouwer 
et  al.  [40]. 


Fig.  6.  The  applied  sinusoidal  voltage  perturbation  with  amplitude  of  0.1  V  under 
different  impedance  frequencies  (10.0-0.01  Hz)  during  one  sinusoidal  cycle. 


140 


M.Y.  Ramandi  et  al.  /  Journal  of  Power  Sources  231  (2013)  134-145 


3.  Numerical  study 

The  computational  domain  is  divided  into  a  number  of  control 
volumes  using  ANSYS  ICEM  CFD  12.0.1.  The  Finite  Volume  based 
commercial  software,  ANSYS  FLUENT  12.0.1,  is  employed  to  dis¬ 
cretize  and  solve  the  set  of  governing  equations.  These  equations 
are  strongly  coupled  through  the  source  terms  resulting  from  the 
electrochemical  reactions.  A  set  of  under-relaxation  techniques  is 
developed  to  handle  the  divergence  difficulties.  Furthermore,  the  C 
programming  language  is  used  to  develop  a  supplemental  code  in 
order  to  include  several  capabilities  (e.g.  non-standard  transport 
equations,  source  terms,  temperature  dependent  properties,  etc.)  to 
ANSYS  FLUENT  12.0.1.  The  SIMPLE  algorithm  is  implemented  for  the 
coupling  between  the  pressure  and  velocity  field.  An  algebraic 
multi-grid  (AMG)  method  with  a  Gauss— Seidel  type  smoother  is 
used  to  accelerate  the  convergence.  A  strict  convergence  criterion 
with  a  residual  of  10-10  is  used  for  all  variables.  The  solution  pro¬ 
cedure  is  illustrated  in  Fig.  2. 

4.  Results  and  discussion 


model  was  employed  to  check  the  validity  of  the  results  by  con¬ 
ducting  a  comparison  with  an  experimental  study  and  satisfactory 
agreement  was  achieved  (Fig.  4). 

4.2.  Non-linear  impedance  load  change 

The  various  transport  phenomena  in  a  molten  carbonate  fuel 
cell  have  a  broad  range  of  time  scales.  It  is  acknowledged  that  the 
electrochemical  reactions  at  the  three-phase  boundaries  in  the 
anode  and  cathode,  heat  conduction  in  solid  electrodes  and  elec¬ 
trolyte,  convective  form  of  heat  and  mass  transfer  in  gas  channels, 
and  charge  transport  processes  at  the  electrodes  and  inside  the 
electrolyte  are  the  dominant  processes  which  take  place  inside  an 
MCFC.  To  be  able  to  elucidate  the  phase  shifts  and  time  scales  for 
different  processes  occurring  in  the  MCFCs,  a  transient  technique  is 
incorporated  into  the  mathematical  model. 

In  this  approach,  a  sinusoidal  (i.e.  non-linear)  voltage  perturba¬ 
tion  is  applied  to  assess  the  resulting  harmonic  response  of  current 
density,  outlet  mass  flow  rates,  heat  transfer  rates,  mass  fractions 
and  temperatures.  The  general  form  of  the  applied  voltage  is 


4.1.  Model  validation 


V(t)  =  Vb+Asin(27r/t) 


(30) 


First  of  all,  for  a  specific  case,  available  experimental  data  was 
used  to  adjust  model  input  parameters  (Fig.  3).  Then,  the  validated 


This  sinusoidal  functional  serves  as  the  boundary  condition  for  the 
solid  phase  potential  at  the  cathode  land  boundaries  (in  Equation  (28)): 


—  Anode  Volumetric  Reaction  Rate  (mol  m'*s) 

•  -  Cathode  Volumetric  Reaction  Rate  (mol  m4s) 
Anode  Outlet  Mass  Flow  Rate  (kg  s  ’) 

—  Cathode  Outlet  Mass  F  low  Rate  (kg  s  ’) 

Sb 


<§> 

&  V 

<o  c$o  w 
y  jc 

c> 

V  TO 

<b'  r  a 

£  | 

•O  o 

g  2 

T3 

<&  © 


O 


Anode  Volumetric  Reaction  Rate  (mol  m 1  s) 


0.04  0.06 

Time  (s) 


M' 

&  ' 
J<b 
1  # 

# 


-4<o  cj 
- 

<b‘ 


4p  o 


(a)  10  Hz 

-  Anode  Volumetric  Reaction  Rate  (mol  m4  s) 

-  Cathode  Volumetric  Reaction  Rate  (mol  m4s) 

-  Anode  Outlet  Mass  F  low  Rate  (kg  s  ') 

-  Cathode  Outlet  Mass  Flow  Rate  (kg  s  ’)  (v. 

_ <v  ~  r^/- 


(b)  1  Hz 


£ 

■e  3  <4- 

&  1  :>  ' 
-  «  JVCV 


<o- ^  8 

Of 

CS5  3 
Oy  3 

o 

| 

ra 

<b*  O 


-  Anode  Volumetric  Reaction  Rate  (mol  m4  s) 

- Cathode  Volumetric  Reaction  Rate  (mol  m4  s) 

- Anode  Outlet  Mass  Flow  Rate  (kg  s  ’) 

- Cathode  Outlet  Mass  Flow  Rate  (kg  s  ’)  ^ 


20  40  60  80  100 

Time  (s) 


£ 

$  S> 

<0  3 

>  & 
* 
i 

*>&  S 

&  s 


<b’  ® 

& 

<c 

o 

v 

<o' 


&  % 


(c)  0.1  Hz 


(d)  0.01  Hz 


Fig.  7.  Dynamic  response  of  the  electrodes  volumetric  reaction  rates  and  gas  flow  rates  at  gas  channel  outlets  corresponding  to  sinusoidal  voltage  perturbation  during  one  si¬ 
nusoidal  cycle  and  over  a  wide  range  of  impedance  frequencies  (10.0-0.01  Hz). 
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=  m 


(31)  /  =  o.OOl,  0.01,  0.1,  1.0,  10.0  Hz 


In  this  equation,  Vb  is  the  base  voltage,  A  is  the  amplitude  of  the 
voltage  perturbation, /is  the  frequency  of  the  sine  wave  and  t  is  the 
operating  time. 

Using  the  boundary  condition  described  in  Equations  (30)  and 
(31)  and  keeping  all  other  input  quantities  constant,  several  sim¬ 
ulations  were  carried  out  in  transient  mode.  The  results  demon¬ 
strate  the  non-linear  dynamic  responses  of  the  unit  cell  variables. 

Generally,  if  the  amplitude  of  the  applied  perturbation  is  small 
enough,  the  dynamic  responses  are  also  sinusoidal  (although  this 
may  not  be  applicable  for  the  energy  (heat)  transport  process  in 
MCFCs).  Accordingly,  the  dynamic  response  of  an  affected  field 
variable,  0,  is 

0  =  05  +A@sin(27r/t  +  6(f))  (32) 

where  0(f)  is  the  phase  difference  between  the  voltage  and  the 
response  function.  Obviously,  for  a  purely  resistive  behaviour  6  is 
zero. 

In  this  study,  the  following  parameters  are  utilized  to  analyze 
the  non-linear  dynamic  responses: 

Vb  =  0.7  V 

A  =  0.1  V 


Fig.  5  demonstrates  the  sinusoidal  perturbation  for  all  fre¬ 
quencies  over  one  periodic  cycle.  In  this  figure,  the  horizontal  axis  is 
the  non-dimensional  time  which  is  defined  as: 


Non  -  dimensional  time 


Time(s) 

One  periodic  cycle(s) 


(33) 


By  this  definition,  the  applied  sinusoidal  voltage  versus  non- 
dimensional  time  follows  identical  shapes  for  all  frequencies. 
Accordingly,  0.50  along  the  horizontal  axis  corresponds  to  500  s, 
50  s,  5  s,  0.5  s  and  0.05  s  for  the  impedance  frequencies  of  0.001  Hz, 
0.01  Hz,  0.1  Hz,  1  Hz  and  10  Hz,  respectively. 

Hereafter,  Fig.  5  along  with  the  resulting  dynamic  responses  will 
be  utilized  to  determine  the  phase  shifts  among  the  different  var¬ 
iables  and  the  applied  perturbation.  Nevertheless,  the  characteristic 
time  scale  of  each  process  will  also  be  assessed. 

The  most  critical  parameter  in  the  dynamic  response  analysis  of 
fuel  cells  is  the  average  cell  current  density.  The  corresponding  time 
scale  characterizes  the  charge  transport  in  both  solid  and  electro¬ 
lyte  phase  and  can  be  evaluated  by  [41]: 


L  charge 
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Fig.  8.  Dynamic  response  of  the  gas  species  mass  fractions  corresponding  to  sinusoidal  voltage  perturbation  during  one  sinusoidal  cycle  and  over  a  wide  range  of  impedance 
frequencies  (10.0-0.01  Hz). 
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where  8  is  the  electrode  thickness  and  C  is  the  capacity  which  ranges 
from  1.5  to  55  F  m”2  [42].  Using  the  model  input  parameters,  the 
time  scale  of  the  charge  transport  process  is  evaluated  to  be  3.6E-06 
to  1.5E-04  s.  These  values  imply  that  the  charge  transport  in  MCFCs  is 
very  fast.  Using  the  evaluated  time  step  size,  the  transient  technique 
is  employed  to  provide  an  improved  understanding  of  the  charge 
transport  process.  Fig.  6  shows  the  dynamic  responses  of  the  unit  cell 
average  current  density  during  one  sinusoidal  cycle  under  four  dif¬ 
ferent  frequencies  (0.01-10  FIz).  It  may  be  observed  that  the 


resulting  dynamic  responses  for  all  impedance  frequencies  follow 
a  sinusoidal  trend  with  no  noticeable  phase  shifts.  This  simply  en¬ 
tails  that  the  dynamic  response  of  the  cell  current  density  is  inde¬ 
pendent  of  the  perturbation  frequency.  Technically,  this  occurs  only 
when  the  associated  time  scale  is  extremely  small. 

The  next  parameters  of  interest  are  the  anode  and  cathode 
volumetric  reaction  rates  which  are  demonstrated  in  Fig.  7  for 
various  frequencies.  This  figure  exhibits  a  pattern  similar  to  the 
average  current  density  with  no  phase  shift.  It  seems  as  if  the 
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Fig.  9.  Dynamic  response  of  electrolyte  and  electrodes’  temperature  and  outlet  energy  corresponding  to  sinusoidal  voltage  perturbation  during  one  sinusoidal  cycle  and  over  a  wide 
range  of  impedance  frequencies  (10.0-0.001  Hz). 
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electrochemical  reaction  rates  respond  to  the  voltage  change 
instantaneously.  This  is  expected  because  the  charge  transport 
process  which  exhibited  an  extremely  small  time  scale  is  coupled 
with  the  electrochemical  reactions.  In  other  words,  the  ionic  and 
electronic  charges  are  generated  and  consumed  on  the  same  time 
scales  as  the  charge  transport  process. 

In  order  to  identify  the  dynamic  behaviour  of  the  mass  transport 
process,  two  types  of  parameters  are  considered.  One  is  the  average 
bulk  value  of  mass  flow  rates  at  gas  channel  outlets  (shown  in 
Fig.  7),  and  the  second  (Fig.  8)  is  the  local  value  of  the  gaseous 


species  mass  fractions  in  the  anode  and  cathode  (at  Points  1  and  3, 
shown  in  Fig.  1 ).  The  approximate  characteristic  time  scale  is  then 
evaluated  as 


T'mass  —  £)eff’ 

assuming  that  diffusion  within  the  electrode  dominates  the  mass 
transport. 


(b) 


(c) 


(d) 


Fig.  10.  Time-extended  dynamic  response  of  average  current  density,  anode  reaction  rate,  AGC  outlet  mass  flow  rate  (at  10  Hz)  and  anode  temperature  (at  1  Hz). 
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By  substituting  the  parameter  values,  the  time  scale  of  the  gas 
transport  process  is  evaluated  to  be  about  0.2-1  s.  This  value  in¬ 
dicates  that  the  gas  transport  process  is  slower  than  the  electro¬ 
chemical  reactions  and  the  charge  transport  process.  Fig.  7  confirms 
this  conclusion  by  showing  obvious  phase  shifts  for  both  AGC  and 
CGC  outlet  mass  flow  rates  at  the  frequency  of  10  Hz.  These  phase 
shifts  then  decrease  and  vanish  at  a  frequency  of  1  Hz.  One  can  justify 
the  fading  of  the  phase  shift  by  taking  into  account  the  sluggishness 
of  the  variation  in  the  operating  condition,  when  the  oscillation 
frequency  is  small.  In  this  situation,  the  unit  cell  has  more  time  to 
reach  its  local  equilibrium.  Therefore,  the  phase  shift  disappears 
when  the  impedance  frequency  diminishes.  Even  though  the  AGC 
and  CGC  outlet  mass  flow  rates  (bulk  variables)  are  reasonable  pa¬ 
rameters  for  analysis  of  the  mass  transport  characteristic  time  scales, 
a  number  of  local  variables  are  also  selected  to  present  an  enhanced 
understanding  of  mass  transport  process.  Results  are  summarized  in 
Fig.  8,  confirming  the  same  time  scales  as  that  of  the  bulk  variables. 

Finally,  the  characteristic  time  scale  of  the  energy  (heat)  trans¬ 
port  process  is  examined.  It  can  be  evaluated  as 


Tenergy  =  (36) 

It  is  worthwhile  to  mention  that  there  are  two  major  energy 
transport  mechanisms  in  an  MCFC.  One  is  the  convective  heat 
transfer  in  gas  channels,  and  the  other  is  the  heat  conduction  in 
electrodes  and  the  electrolyte.  The  overall  energy  transport  process  is 
controlled  by  the  slowest  mechanism  which  has  the  largest  time 
scale.  Since  the  heat  capacity  of  the  solid  material  is  much  higher  than 
the  heat  capacity  of  the  gas  mixture,  the  parameters  in  Equation  (36) 
are  replaced  by  the  solid  material  properties.  Thus,  the  characteristic 
time  scale  is  found  to  be  of  the  order  of  1000  s.  This  value  shows  that 
the  energy  transport  process  is  three  orders  of  magnitude  slower 
than  the  mass  transport  process.  This  notable  conclusion  facilitates 
the  process  of  choosing  the  appropriate  time  step  size.  Transient 
simulations  must  be  started  with  an  extremely  small  time  step  size 
(1.0E-6.0  s)  which  can  be  increased  up  to  1  s  when  the  operating  time 
passes  the  mass  transport  time  scale.  This  way,  solution  accuracy  is 
secured  along  with  optimal  computational  expense. 

Similar  to  the  mass  transport  dynamic  response,  both  local  and 
bulk  variables  are  chosen  to  analyze  the  phase  shifts  for  energy 
transport  characteristics.  Fig.  9  shows  the  resulting  dynamic 
response  of  the  local  temperature  (at  Points  1,  2  and  3  in  Fig.  1 )  and 
the  average  bulk  value  of  total  energy  at  gas  channel  outlets.  It  may 
be  observed  that  there  is  more  complexity  in  the  dynamic  responses 
corresponding  to  the  energy  transport  process.  First  of  all,  the 
electrolyte  temperature  does  not  exhibit  a  simple  sinusoidal  shape 
at  high  frequencies  (10  Hz  and  1  Hz  in  Fig.  9a  and  b).  This  underlines 
that  the  heat  conduction  process  in  the  electrolyte  has  the  largest 
time  scale  and  hence  is  the  slowest  thermal  process  in  the  unit  cell. 
One  may  associate  this  with  the  heat  conductivity  or  the  heat  ca¬ 
pacity  of  the  electrolyte.  However,  the  sluggishness  of  this  process 
can  only  be  justified  by  Equation  (36)  which  shows  the  largest  value 
for  the  electrolyte  time  scale.  Nonetheless,  as  the  impedance  fre¬ 
quency  decreases,  the  electrolyte  temperature  dynamic  response 
recovers  the  sinusoidal  shape.  Additionally,  there  are  obvious  phase 
shifts  for  all  frequencies  (Fig.  9a— d)  regarding  the  dynamic  re¬ 
sponses  of  the  four  variables.  Therefore,  another  simulation  was 
carried  out  to  examine  the  dynamic  response  with  a  much  smaller 
frequency  (0.001  Hz).  Results  are  shown  in  Fig.  9e.  Even  at  this  small 
frequency,  the  thermal  characteristics  of  the  unit  cell  do  not  show 
a  rapid  response  to  the  voltage  perturbation.  Also,  the  local  and  bulk 
values  of  the  investigated  variables  and  parameters  demonstrate 
different  phase  shifts  indicating  their  different  characteristic  time 
scales.  The  phase  shifts  can  be  predicted  to  disappear  at  frequencies 


Fig.  11.  Linear  dynamic  response  of  the  anode  temperature  to  a  small  change  in  the 
operating  voltage. 


under  0.0005  Hz.  This  conclusion  reveals  the  large  time  scale  of  the 
energy  transport  process  to  be  2500  s.  This  great  hindrance  of  the 
dynamic  response,  in  comparison  with  the  mass  and  charge 
transport  processes,  is  primarily  because  of  the  high  thermal  ca¬ 
pacity  of  the  fuel  cell.  This  feature  can  be  beneficial  for  the  design  of 
control  systems  for  MCFC-hybrid  plants. 

Figs.  7-9  show  that  the  majority  of  the  chemical  and  physical 
processes  in  an  MCFC  occur  in  the  first  second  of  operation.  Only 
heat  transport  takes  longer  to  reach  its  local  equilibrium.  Therefore, 
the  same  impedance  simulation  was  extended  for  20  cycles  at 
a  frequency  of  10  Hz  to  study  the  dynamic  responses  of  the  fast  and 
very  fast  processes.  This  large  frequency  is  required  to  capture  the 
transient  period.  However,  for  the  sake  of  comparison,  only  a  few 
variables  are  selected  for  the  extended  time  study.  Results  are 
summarized  in  Fig.  lOa-c.  It  may  be  observed  that  the  average  cell 
current  density  (Fig.  10a)  and  anode  reaction  rate  (Fig.  10b)  respond 
to  the  voltage  fluctuation  instantaneously  without  any  obvious 
phase  shift.  On  the  other  hand,  the  response  of  the  AGC  outlet  mass 
flow  rate  (Fig.  10c)  is  relatively  slower.  A  0.2  s  transition  time  is 
required  for  the  outlet  flow  rate  to  reach  the  steady-state  condition. 
The  last  part  of  Fig.  10  (lOd)  demonstrates  the  anode  temperature 
response  for  the  first  50  cycles  at  an  impedance  frequency  of  1  Hz. 
This  small  frequency  is  necessary  because  the  energy  transport 
process  is  a  much  slower  process.  Fig.  lOd  shows  that  the  oscillation 
amplitude  of  the  anode  temperature  increases  instantly,  and  has 
a  continuous  growing  trend  throughout  the  operation  time.  This 
simulation  was  further  continued  for  the  next  1000  s  and  the  same 
rising  trend  was  observed.  Due  to  the  enormous  number  of  cycles 
over  the  extended  operation  time,  the  resulting  graph  is  not  read¬ 
able  and  helpful.  Thus,  it  is  not  included  in  this  study. 

In  its  place,  a  new  simulation,  with  a  linear  voltage  change 
during  the  initial  phase  of  the  operation,  was  carried  out  and  the 
results  are  demonstrated  in  Fig.  11.  This  figure  shows  the  anode 
temperature  variation  over  a  time  period  of  4500  s.  It  confirms  the 
previous  conclusion  (from  Fig.  9e)  that  the  transition  period  of  the 
energy  transport  process  is  about  2500  s. 

5.  Conclusions 

A  three-dimensional  transient  mathematical  model  for  a  MCFC 
was  developed  to  study  the  dynamic  response  of  the  fuel  cell. 
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ANSYS  FLUENT  12.0.1,  a  finite  volume  based  commercial  software 
package,  was  utilized  to  solve  the  system  of  partial  differential 
equations.  A  sinusoidal  impedance  approach  was  employed  to 
identify  the  characteristic  time  scales  for  the  major  dynamic 
transport  processes.  Thus,  the  electrochemical  reactions,  charge 
transport,  mass  transport  and  energy  transport  processes  were 
analyzed  to  find  the  cell  response  during  a  sinusoidal  voltage 
change  and  over  a  wide  range  of  impedance  frequencies.  The  cor¬ 
responding  time  scales  were  verified  numerically  with  distinctive 
order  of  magnitudes,  giving  us  confidence  in  our  model  and  nu¬ 
merical  solutions.  The  majority  of  the  physical  and  electrochemical 
processes  occur  in  the  first  second  of  the  operation.  The  anodic  and 
cathodic  electrochemical  reactions  along  with  the  charge  transport 
process  were  found  to  be  the  fastest  processes  in  the  unit  cell  with 
time  scales  about  10  6-10  4  s.  The  time  scale  of  the  gas  transport 
process  was  evaluated  to  be  in  order  of  1  s,  while  the  energy 
transport  process  exhibited  a  time  scale  larger  than  1000  s.  This  is 
in  agreement  with  experimental  data  provided  internally  by 
Enbridge  Inc.  [49].  Hence,  a  variable  step  size  algorithm  was 
employed  to  increase  the  efficiency  of  our  computations. 

Having  established  a  functional  transient  model  and  code, 
future  work  will  employ  the  developed  model  along  with  the  dy¬ 
namic  analysis  result  to  design  more  robust  process  control  strat¬ 
egies  for  MCFC-hybrid  systems.  This  work  will  also  be  used  as  a  first 
step  to  design  a  diagnosis  tool  in  thermal-stress  analysis  to  deter¬ 
mine  degradation  effects,  and  may  help  to  extend  the  life  cycle  of 
the  fuel  cell. 
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